function [tvhat, est_ols, residuals, test_ols, mdl] = compute_ols_inference(...
    y, x, xi, nw_lag_trunc, test_beta, asy_cov_matrix)
    % compute_ols_inference computes the OLS estimator, residuals, 
    % t-statistics, and standard errors for the OLS coefficients
    % normalized by a given normalizing matrix and corrected for 
    % heteroskedasticity and autocorrelation, considering 
    % the deterministic time trend in the independent variable.
    %
    % Inputs:
    % y              - t x 1 column vector of dependent variables
    % x              - t x k matrix of independent variables (should include the intercept)
    % xi             - scalar parameter that affects the structure of the covariance matrix
    % nw_lag_trunc   - scalar, truncation parameter for the Newey-West kernel estimator
    % test_beta      - k x l matrix of the null hypothesis values for the OLS coefficients (for up-to l sets of t-test)
    % asy_cov_matrix - boolean, whether to use the asymptotic covariance matrix or finite sample version
    %
    % Outputs:
    % tvhat         - k x 1 vector of normalized standard errors
    % est_ols       - k x 1 vector of OLS coefficient estimates
    % residuals     - t x 1 vector of residuals from the OLS model
    % test_ols      - k x l vector of t-statistics for the OLS coefficients based on test_beta
    % mdl           

    % number of observations (time periods)
    t = size(y, 1);

    % number of independent variables
    k = size(x, 2);

    % Ensure test_beta is the same dimension as est_ols
    if size(test_beta, 1) ~= k
        error('Dimension mismatch between est_ols and test_beta');
    end

    % Select the appropriate covariance matrix based on the input flag `asy_cov_matrix`
    if asy_cov_matrix
        % Asymptotic covariance matrix
        Xi_11 = [1, 1/2; 1/2, 1/3];
        Xi_12 = [sqrt(1-xi), sqrt(1-xi)/2; sqrt(1-xi)*(1+xi)/2, sqrt(1-xi)*(2+xi)/6];
        Xi_21 = Xi_12';
        Xi_22 = Xi_11;
    else
        % Finite sample covariance matrix
        Xi_11 = [1, (t+1)/(2*t); (t+1)/(2*t), ((2*t+1)*(t+1))/(6*t^2)];
        Xi_12 = [sqrt(1-xi), (1+(1-xi)*t)/(2*t*(1-xi)^(1/2)); ...
                 ((1-xi)^(1/2)*((1+xi)*t+1))/(2*t), ...
                 (1+(1-xi)*(1+xi)*t^2+3*t)/(6*t^2*(1-xi)^(1/2))];
        Xi_21 = Xi_12';
        Xi_22 = [1, ((1-xi)*t+1)/(2*(1-xi)*t); ((1-xi)*t+1)/(2*(1-xi)*t), ...
                 (2*(1-xi)^2*t^2+3*(1-xi)*t+1)/(6*t^2*(1-xi)^2)];
    end

    % Construct the full covariance matrix
    cov_mat = [Xi_11, Xi_12; Xi_21, Xi_22];

    % Calculate the inverse of the covariance matrix
    inv_cov_mat = inv(cov_mat);

    % Normalize the covariance matrix using the appropriate normalizing matrix
    norm_mat = [sqrt(t), t^(3/2), sqrt(t*(1-xi)), (t*(1-xi))^(3/2)];

    % Estimate the OLS model and retrieve coefficients (without intercept in x)
    mdl = fitlm(x, y, 'intercept', false);
    est_ols = reshape(mdl.Coefficients.Estimate, [k, 1]);

    % Compute residuals from the model
    residuals = mdl.Residuals.Raw;

    % Calculate the long-run variance using the Newey-West kernel estimator
    long_run_variance = kernelEstimator(ones(size(residuals)), residuals, ...
        nw_lag_trunc, 'NW') * t;

    % Compute the covariance matrix of the coefficients
    tvb = long_run_variance .* inv_cov_mat;

    % Compute normalized standard errors
    tvhat = sqrt(diag(tvb)) ./ norm_mat';

    % Compute the t-statistics for OLS coefficients based on the null hypothesis (test_beta)
    test_ols = (est_ols - test_beta) ./ tvhat;

end